ENMTools

R build status

This package implements various tests, visualizations, and metrics for use with environmental niche models (ENMs) and species distribution models (SDMs).

Citation

Warren, D.L., N. Matzke, M. Cardillo, J. Baumgartner, L. Beaumont, N. Huron, M. Simões, Teresa L. Iglesias, and R. Dinnage. 2019. ENMTools (Software Package). URL: https://github.com/danlwarren/ENMTools. doi:10.5281/zenodo.3268814

Installation

At present, ENMTools is downloadable from https://github.com/danlwarren/ENMTools. There are multiple ways to download it. The easiest is to use devtools and install from GitHub.

Installing from GitHub using devtools

Run the following code from your R console:

install.packages("devtools")
library(devtools)
install_github("danlwarren/ENMTools")
library(ENMTools)

Install from zip file

A zipped version of the package is available at https://github.com/danlwarren/ENMTools/archive/master.zip. To install from the zip file, download a copy of it to your system. Once it’s finished downloading, type the following (where PATH is the path to the zip file):

install.packages("devtools")
library(devtools)
install_local("PATH")
library(ENMTools)

Installing extras

ENMTools uses functionality from a LOT of other R packages, and it’s possible that you don’t want to install them all. For that reason many of the packages are not automatically installed with ENMTools, but instead “suggested”. If you want to install all of the suggested packages, we have a function for that. You should only need to use it after you first install ENMTools or update R. If you choose not to install the extra packages, you will get warnings when you try to use functions that require them.

install.extras()

Interacting with ENMTools

Creating enmtools.species objects

We’re going to load in some environmental data. You can do this from local rasters, like so:

env.files <- list.files(path = "./env_pca/", pattern = "pc", full.names = TRUE)
env <- stack(env.files)
names(env) <- c("pc1", "pc2", "pc3", "pc4")
env <- setMinMax(env)

Or you can load them from the internet using the raster package’s getData() function.

library(raster)
env <- raster::getData('worldclim', var='bio', res=10)
env <- crop(env, extent(-10, 17, 39, 48))
plot(env[[1]])

ENMTools is primarily designed to examine patterns of similarity and difference between ENMs for different species. In order to simplify interactions with the functions in ENMTools, you need to put your data for each of your species into an enmtools.species object. You can create and view an empty enmtools.species object just by typing:

monticola <- enmtools.species()
monticola
## 
## 
## Range raster not defined.
## 
## Presence points not defined.
## 
## Background points not defined.
## 
## Species name not defined.

You can add bits of it when the object is created:

monticola.path <- paste(system.file(package="ENMTools"), "/monticola.csv", sep='')
monticola <- enmtools.species(species.name = "monticola", 
                            presence.points = read.csv(monticola.path))
monticola$range <- background.raster.buffer(monticola$presence.points, 50000, mask = env)
monticola$background.points <- background.points.buffer(points = monticola$presence.points,
                                                   radius = 20000, n = 1000, mask = env[[1]])

Or you can add data to this object after it’s created:

names(monticola)
monticola$species.name <- "monticola"
monticola$presence.points <- read.csv(monticola.path)
monticola$range <- background.raster.buffer(monticola$presence.points, 50000, mask = env)
monticola$background.points <- background.points.buffer(points = monticola$presence.points,
                                                   radius = 20000, n = 1000, mask = env[[1]])

It’s always a very good idea to run the check.species() function on an enmtools.species object after you build or modify it.

monticola <- check.species(monticola)

And now we can take a look at our species!

interactive.plot.enmtools.species(monticola)

However, ENMTools also contains some sample data. It contains an enmtools.clade object called “iberolacerta.clade”, which holds several enmtools.species objects including an already-built version of Iberolacerta monticola. It also has some low-resolution Worldclim data that we can use to demonstrate functions. We’ll pull two of those species out now.

data(iberolacerta.clade)
data(euro.worldclim)
monticola <- iberolacerta.clade$species$monticola
cyreni <- iberolacerta.clade$species$cyreni
env <- euro.worldclim

Building an ENM

ENMTools contains functions to simplify the ENM construction process. Using enmtools.species objects and the correct modeling commands, we can build models very quickly. These commands are primarily wrappers to dismo model construction and projection functions, and at present are only available for GLM, Maxent, Domain, and Bioclim models. One of the nice bits about this setup is that it allows enmtools to automatically generate suitability maps, do model evaluation, and plot the marginal suitability of habitat for each variable separately.

Before we get started modeling, however, we might want to trim our predictor set so as to reduce collinearity. ENMTools contains a couple of handy functions for this.

raster.cor.matrix(env)
##              bio1         bio2        bio3        bio4       bio5        bio6
## bio1   1.00000000  0.208715969  0.36766171 -0.21357041  0.8853513  0.93295069
## bio2   0.20871597  1.000000000  0.58362743  0.30918934  0.5417192 -0.01723188
## bio3   0.36766171  0.583627431  1.00000000 -0.56601419  0.3052211  0.42521254
## bio4  -0.21357041  0.309189341 -0.56601419  1.00000000  0.1754180 -0.52131150
## bio5   0.88535133  0.541719174  0.30522110  0.17541802  1.0000000  0.68751772
## bio6   0.93295069 -0.017231880  0.42521254 -0.52131150  0.6875177  1.00000000
## bio7  -0.05329392  0.709020870 -0.14900115  0.87998888  0.4016817 -0.38884596
## bio8   0.24945265  0.068562982 -0.16262988  0.30400935  0.2862188  0.09285337
## bio9   0.73849291  0.183148800  0.42195497 -0.37866158  0.6400464  0.77343845
## bio10  0.96045067  0.311324869  0.21975162  0.06247355  0.9606340  0.80726704
## bio11  0.95921479  0.101719882  0.48107648 -0.47305206  0.7533941  0.98896057
## bio12 -0.60846674 -0.482731496 -0.31488066 -0.07961919 -0.7137513 -0.48704454
## bio13 -0.41433559 -0.510317568 -0.25626002 -0.19124801 -0.5722885 -0.28138429
## bio14 -0.72871981 -0.338229901 -0.37618609  0.16712793 -0.7297745 -0.67564389
## bio15  0.45421260 -0.007931916  0.11838593 -0.21759678  0.3821863  0.47445069
## bio16 -0.44597594 -0.506032477 -0.26995545 -0.16924912 -0.5920746 -0.31485467
## bio17 -0.70338945 -0.343806007 -0.33935275  0.10832583 -0.7243068 -0.63769389
## bio18 -0.81091065 -0.331213024 -0.47078459  0.28348226 -0.7776879 -0.80126262
## bio19 -0.12414066 -0.396445432  0.07348362 -0.48121938 -0.3430162  0.07715578
##              bio7        bio8       bio9       bio10       bio11       bio12
## bio1  -0.05329392  0.24945265  0.7384929  0.96045067  0.95921479 -0.60846674
## bio2   0.70902087  0.06856298  0.1831488  0.31132487  0.10171988 -0.48273150
## bio3  -0.14900115 -0.16262988  0.4219550  0.21975162  0.48107648 -0.31488066
## bio4   0.87998888  0.30400935 -0.3786616  0.06247355 -0.47305206 -0.07961919
## bio5   0.40168174  0.28621877  0.6400464  0.96063401  0.75339408 -0.71375135
## bio6  -0.38884596  0.09285337  0.7734384  0.80726704  0.98896057 -0.48704454
## bio7   1.00000000  0.24603268 -0.1633540  0.20071994 -0.29134540 -0.29133189
## bio8   0.24603268  1.00000000 -0.2200483  0.33135353  0.12985480 -0.30622476
## bio9  -0.16335400 -0.22004834  1.0000000  0.66175815  0.78861884 -0.46297347
## bio10  0.20071994  0.33135353  0.6617581  1.00000000  0.84907884 -0.65829840
## bio11 -0.29134540  0.12985480  0.7886188  0.84907884  1.00000000 -0.54163405
## bio12 -0.29133189 -0.30622476 -0.4629735 -0.65829840 -0.54163405  1.00000000
## bio13 -0.37121605 -0.19351156 -0.3074562 -0.48860809 -0.32972845  0.90648180
## bio14 -0.07381576 -0.25699492 -0.6265099 -0.71411350 -0.72423144  0.77619752
## bio15 -0.11344874  0.17648853  0.3970693  0.41693023  0.49288586 -0.24313910
## bio16 -0.35410914 -0.21747973 -0.3274923 -0.51511693 -0.36513899  0.91789736
## bio17 -0.11473790 -0.30488093 -0.5552129 -0.70368996 -0.68356379  0.80791339
## bio18  0.02381502  0.03291308 -0.8163353 -0.76592247 -0.82992875  0.76402046
## bio19 -0.53249335 -0.51356967  0.1072073 -0.26914814  0.01511498  0.74717071
##            bio13       bio14        bio15      bio16      bio17       bio18
## bio1  -0.4143356 -0.72871981  0.454212596 -0.4459759 -0.7033895 -0.81091065
## bio2  -0.5103176 -0.33822990 -0.007931916 -0.5060325 -0.3438060 -0.33121302
## bio3  -0.2562600 -0.37618609  0.118385927 -0.2699555 -0.3393528 -0.47078459
## bio4  -0.1912480  0.16712793 -0.217596785 -0.1692491  0.1083258  0.28348226
## bio5  -0.5722885 -0.72977455  0.382186338 -0.5920746 -0.7243068 -0.77768793
## bio6  -0.2813843 -0.67564389  0.474450688 -0.3148547 -0.6376939 -0.80126262
## bio7  -0.3712161 -0.07381576 -0.113448735 -0.3541091 -0.1147379  0.02381502
## bio8  -0.1935116 -0.25699492  0.176488534 -0.2174797 -0.3048809  0.03291308
## bio9  -0.3074562 -0.62650986  0.397069264 -0.3274923 -0.5552129 -0.81633533
## bio10 -0.4886081 -0.71411350  0.416930231 -0.5151169 -0.7036900 -0.76592247
## bio11 -0.3297284 -0.72423144  0.492885864 -0.3651390 -0.6835638 -0.82992875
## bio12  0.9064818  0.77619752 -0.243139096  0.9178974  0.8079134  0.76402046
## bio13  1.0000000  0.48404051  0.130409758  0.9902335  0.5245850  0.58835234
## bio14  0.4840405  1.00000000 -0.743674332  0.5012004  0.9859234  0.84673991
## bio15  0.1304098 -0.74367433  1.000000000  0.1244492 -0.7294802 -0.46226779
## bio16  0.9902335  0.50120043  0.124449153  1.0000000  0.5362200  0.61586400
## bio17  0.5245850  0.98592339 -0.729480228  0.5362200  1.0000000  0.81764699
## bio18  0.5883523  0.84673991 -0.462267795  0.6158640  0.8176470  1.00000000
## bio19  0.7906993  0.31690226  0.128299899  0.7920907  0.3839700  0.17654879
##             bio19
## bio1  -0.12414066
## bio2  -0.39644543
## bio3   0.07348362
## bio4  -0.48121938
## bio5  -0.34301616
## bio6   0.07715578
## bio7  -0.53249335
## bio8  -0.51356967
## bio9   0.10720732
## bio10 -0.26914814
## bio11  0.01511498
## bio12  0.74717071
## bio13  0.79069926
## bio14  0.31690226
## bio15  0.12829990
## bio16  0.79209065
## bio17  0.38396995
## bio18  0.17654879
## bio19  1.00000000

That’s great, but it’s a bit hard to pick variables this way. Let’s try it visually instead.

raster.cor.plot(env)
## $cor.mds.plot

## 
## $cor.heatmap

The raster.cor.plot function gives us two visualizations. One heatmap that colors pairs of predictors by their Pearson correlation coefficient, and one cluster plot that does mds scaling of the predictor variables and then plots them in a two dimensional space so that more correlated predictors are closer to each other. We’re going to make an arbitrary decision to just use three predictors, and to keep those predictors relatively uncorrelated we’ll select predictors that are far apart in this mds plot. Here we’ll choose bio1, bio12, and bio7.

env <- env[[c("bio1", "bio12", "bio7")]]
plot(env)

raster.cor.matrix(env)
##              bio1      bio12        bio7
## bio1   1.00000000 -0.6084667 -0.05329392
## bio12 -0.60846674  1.0000000 -0.29133189
## bio7  -0.05329392 -0.2913319  1.00000000

GLM

GLMs usually require the user to supply a formula, an enmtools.species object, and some environmental data. If your formula is a strictly additive function of all of the environmental layers in env, though, enmtools.glm will build a formula automatically.

monticola.glm <- enmtools.glm(species = monticola, env = env, f = pres ~ bio1 + bio12 + bio7, test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
monticola.glm
## 
## 
## Formula:  presence ~ bio1 + bio12 + bio7
## <environment: 0x7fc19844e278>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   bio1   bio12   bio7   presence
## ---  ----------  ---------  -----  ------  -----  ---------
## 2     -6.036635   43.02531     76    1012    246          1
## 3     -7.679727   40.38852    137    1143    247          1
## 4     -7.790437   40.30959    129    1231    242          1
## 5     -7.473340   43.78935    140     931    179          1
## 6     -6.575039   42.91070     84    1012    247          1
## 7     -5.132756   43.49572    133     822    190          1
## 8     -7.787378   40.39362    137    1143    247          1
## 9     -4.941888   43.35310    128     843    194          1
## 11    -7.645674   40.36543    101    1514    229          1
## 13    -6.990000   42.57000    107     893    253          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5825  -0.7952  -0.5213   0.8570   2.1071  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  9.1832804  2.2491451   4.083 4.45e-05 ***
## bio1        -0.0346144  0.0064428  -5.373 7.76e-08 ***
## bio12       -0.0004071  0.0007446  -0.547    0.585    
## bio7        -0.0211369  0.0050638  -4.174 2.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 576.70  on 719  degrees of freedom
## Residual deviance: 518.09  on 716  degrees of freedom
## AIC: 256.31
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 512 
## AUC            : 0.7242713 
## cor            : 0.2978555 
## max TPR+TNR at : 0.05330141 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 10000 
## AUC            : 0.4187971 
## cor            : -0.020402 
## max TPR+TNR at : 0.402039 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 512 
## AUC            : 0.7627704 
## cor            : 0.2177419 
## max TPR+TNR at : 0.06321427 
## 
## 
## Environment space model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 10000 
## AUC            : 0.4384788 
## cor            : -0.003889003 
## max TPR+TNR at : 0.4272269 
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 0.0318854, 0.9953804  (min, max)
## 
## 
## 
## Notes:

Notice this produces the same formula as:

monticola.glm <- enmtools.glm(species = monticola, env = env, test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
monticola.glm
## 
## 
## Formula:  presence ~ bio1 + bio12 + bio7
## <environment: 0x7fc196136858>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   bio1   bio12   bio7   presence
## ---  ----------  ---------  -----  ------  -----  ---------
## 2     -6.036635   43.02531     76    1012    246          1
## 3     -7.679727   40.38852    137    1143    247          1
## 4     -7.790437   40.30959    129    1231    242          1
## 6     -6.575039   42.91070     84    1012    247          1
## 8     -7.787378   40.39362    137    1143    247          1
## 11    -7.645674   40.36543    101    1514    229          1
## 12    -7.642539   40.36317    101    1514    229          1
## 13    -6.990000   42.57000    107     893    253          1
## 14    -7.950000   43.30000    120    1200    194          1
## 15    -7.100000   42.93000    115     935    234          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6333  -0.7957  -0.5017   0.8428   2.2371  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  8.5864023  2.2279823   3.854 0.000116 ***
## bio1        -0.0345127  0.0064959  -5.313 1.08e-07 ***
## bio12        0.0001694  0.0007326   0.231 0.817144    
## bio7        -0.0209675  0.0050438  -4.157 3.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 576.70  on 719  degrees of freedom
## Residual deviance: 508.74  on 716  degrees of freedom
## AIC: 250.07
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 512 
## AUC            : 0.7367037 
## cor            : 0.3182308 
## max TPR+TNR at : 0.09138636 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 10000 
## AUC            : 0.4265486 
## cor            : -0.01570869 
## max TPR+TNR at : 0.4964714 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 512 
## AUC            : 0.7059796 
## cor            : 0.1702158 
## max TPR+TNR at : 0.00553722 
## 
## 
## Environment space model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 10000 
## AUC            : 0.4062423 
## cor            : -0.01320871 
## max TPR+TNR at : 0.3522604 
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 0.02485564, 0.9967967  (min, max)
## 
## 
## 
## Notes:

If you want a more complicated formula, though (e.g., with interactions or polynomial effects), you’ll need to supply that manually.

monticola.glm <- enmtools.glm(species = monticola, env = env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) * poly(bio12, 2), test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
monticola.glm
## 
## 
## Formula:  presence ~ poly(bio1, 2) + poly(bio7, 2) + poly(bio12, 2) + poly(bio7, 
##     2):poly(bio12, 2)
## <environment: 0x7fc194a981b0>
## 
## 
## Data table (top ten lines): 
## 
##       Longitude   Latitude   bio1   bio12   bio7   presence
## ---  ----------  ---------  -----  ------  -----  ---------
## 1     -5.171215   43.06957     78     917    249          1
## 2     -6.036635   43.02531     76    1012    246          1
## 4     -7.790437   40.30959    129    1231    242          1
## 5     -7.473340   43.78935    140     931    179          1
## 6     -6.575039   42.91070     84    1012    247          1
## 7     -5.132756   43.49572    133     822    190          1
## 10    -7.621731   40.34170    101    1514    229          1
## 11    -7.645674   40.36543    101    1514    229          1
## 12    -7.642539   40.36317    101    1514    229          1
## 13    -6.990000   42.57000    107     893    253          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3555  -0.7288  -0.3902   0.6606   2.7248  
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -0.7349     0.3950  -1.861   0.0628 .  
## poly(bio1, 2)1                  -35.3240     6.2353  -5.665 1.47e-08 ***
## poly(bio1, 2)2                  -34.4578     6.6670  -5.168 2.36e-07 ***
## poly(bio7, 2)1                  -31.1760    12.9027  -2.416   0.0157 *  
## poly(bio7, 2)2                   14.0819    10.3603   1.359   0.1741    
## poly(bio12, 2)1                  19.4431    12.1623   1.599   0.1099    
## poly(bio12, 2)2                  -4.2766     9.3814  -0.456   0.6485    
## poly(bio7, 2)1:poly(bio12, 2)1  114.9801   419.8595   0.274   0.7842    
## poly(bio7, 2)2:poly(bio12, 2)1 -324.6524   278.6833  -1.165   0.2440    
## poly(bio7, 2)1:poly(bio12, 2)2  329.5345   269.2964   1.224   0.2211    
## poly(bio7, 2)2:poly(bio12, 2)2  236.4078   162.8632   1.452   0.1466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 576.70  on 719  degrees of freedom
## Residual deviance: 454.57  on 709  degrees of freedom
## AIC: 240.46
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 512 
## AUC            : 0.7864427 
## cor            : 0.4027581 
## max TPR+TNR at : 0.110028 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 208 
## n absences     : 10000 
## AUC            : 0.6633053 
## cor            : 0.08390611 
## max TPR+TNR at : 0.3260363 
## 
## 
## Proportion of data wittheld for model fitting:  0.2
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 512 
## AUC            : 0.8040114 
## cor            : 0.2623441 
## max TPR+TNR at : 0.2500096 
## 
## 
## Environment space model fit (test data):  class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 10000 
## AUC            : 0.6672923 
## cor            : 0.04552655 
## max TPR+TNR at : 0.5621803 
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 2.220446e-16, 0.9751367  (min, max)
## 
## 
## 
## Notes:

To check out the marginal response functions, you only need to type

monticola.glm$response.plots
## $bio1

## 
## $bio7

## 
## $bio12

These plots present a smoothed estimate of the frequency of different levels of the environemntal variable in the presence data and the background points, along with the estimated relationship between that environmental predictor and the suitability of habitat from the model.

You can also visualize your models and data in a 2D environment space using any pair of layers from your environment stack. These plots hold all non-plotted variables (bio7 in this case) constant at their mean value across all presence points, then vary the plotted variables between the minimum and maximum values in env.

The suit.plot shows you suitability in environment space as a function of your two variables, with brighter colors representing variable combinations predicted to be more suitable. The points represent the occurrence points for your species in that environment space.

The colored raster of the background.plot shows you the density of background points in environment space, while the white points again represent your occurrence points in environment space.

visualize.enm(monticola.glm, env, layers = c("bio1", "bio12"), plot.test.data = TRUE)
## $background.plot
## Warning: Removed 396 rows containing missing values (geom_raster).

## 
## $suit.plot

GAM, Bioclim, Domain, and Maxent

The procedure for building Bioclim, Domain, and Maxent models is similar to the procedure for GLMs, with the exception that you do not need to pass a formula to the model function for Maxent, Domain, and Bioclim models. Note that running Maxent models requires a bit of extra setup; see dismo documentation for details.

monticola.gam <- enmtools.gam(monticola, env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) * poly(bio12, 2), test.prop = 0.2)
monticola.dm <- enmtools.dm(monticola, env, test.prop = 0.2)
monticola.bc <- enmtools.bc(monticola, env, test.prop = 0.2)
monticola.mx <- enmtools.maxent(monticola, env, test.prop = 0.2)

Metrics: breadth, correlation, and overlap

ENMTools provides a number of metrics for ENMs and for similarities between ENMs. These include measures of niche breadth, based on Levins(1968). An important caveat when interpreting these metrics is that they are driven to some (variable) extent by the availability of different combinations of environmental predictors. As such they are more accurately interpreted as a measurment of the smoothness of the geographic distribution of suitability scores than as an estimate of the breadth of the fundamental niche; an organism with narrow fundamental niche breadth that nonetheless encompasses a set of environmental conditions that is quite common will have a high breadth when measured using ENMs, while having a low breadth in environment space.

raster.breadth(monticola.glm)
## $B1
## [1] 0.9448514
## 
## $B2
## [1] 0.5414817

ENMTools also provides metrics for measuring similarity between ENMs. These include Schoener’s D (Schoener 1968), I (Warren et al. 2008), and the Spearman rank correlation coefficient between two rasters. While D and I are commonly used in the ENM literature, they may tend to overestimate similarity between ENMs when many grid cells are of similar values (e.g., when two species prefer different habitat but the region contains a great deal of habitat that is unsuitable for both).

monticola.glm <- enmtools.glm(species = monticola, env = env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) + poly(bio12, 2), test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
cyreni.glm <- enmtools.glm(species = cyreni, env = env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) + poly(bio12, 2), test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
raster.overlap(monticola.glm, cyreni.glm)
## $D
## [1] 0.6836264
## 
## $I
## [1] 0.8944661
## 
## $rank.cor
## [1] 0.2951514

A new feature of the R version of ENMTools is that you can now use these same metrics in the n-dimensional space of all combinations of environmental variables, instead of restricting your measures of model similarity to those sets of conditions that appear in the training region. This is done by repeatedly drawing Latin hypercube samples from the space of all possible combinations of environmental variables given the min and max of each variable within the training region. ENMTools continues to draw samples until subsequent iterations differ by less than a specified tolerance value. Lower tolerance values result in more precise estimates of overlap, but can take much longer to calculate.

monticola.glm <- enmtools.glm(species = monticola, env = env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) + poly(bio12, 2), test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
cyreni.glm <- enmtools.glm(species = monticola, env = env, f = pres ~ poly(bio1, 2) + poly(bio7, 2) + poly(bio12, 2), test.prop = 0.2)
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
env.overlap(monticola.glm, cyreni.glm, env, tolerance = .001)
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## $env.D
## [1] 0.9206785
## 
## $env.I
## [1] 0.9894104
## 
## $env.cor
## [1] 0.9886251
## 
## $env.D.plot

## 
## $env.I.plot

## 
## $env.cor.plot

The plots that come out of these environment space functions are used for diagnosing convergence of the overlap/breadth metric. Ideally what you want is a relationship between the metric and the number of samples that shows no clear directional trend.

Hypothesis testing

Niche identity or equivalency test

In this example, we will run a niche identity (also called equivalency) test, as in Warren et al. 2008. This test takes the presence points for a pair of species and randomly reassigns them to each species, then builds ENMs for these randomized occurrences. By doing this many times, we can estimate the probability distribution for ENM overlap between species under the null hypothesis that the two species’ occurrences in the environment are effectively a random draw from the same underlying distribution. Note that niche evolution is only one of many reasons why two species’ realized environmental distributions might cause departures from this null hypothesis. See Warren et al. 2014 for details.

To run an identity test, we need to decide what type of models we will build, how many replicates we will run, and (in the case of GLM and GAM) a model formula to use for empirical models and the Monte Carlo replicates. The resulting object contains the replicate models, p values, and plots of the results. Typically idenity tests are run with at least 99 replicates, but we are using a smaller number here for the sake of execution time.

NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the identity test are conducted with pseudoabsence points pooled for the two species being compared.

id.glm <- identity.test(species.1 = monticola, species.2 = cyreni, env = env, type = "glm", nreps = 4)
id.glm
## 
## 
##  
## 
## Identity test monticola vs. cyreni
## 
## objectentity test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I    rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  ----------  ----------  ----------  ----------
## empirical    0.3244487   0.5929443   0.1049357   0.3668837   0.6148383   0.1351514
## rep 1        0.9718264   0.9992093   0.9883338   0.9543520   0.9973478   0.9879216
## rep 2        0.9380929   0.9961221   0.9792829   0.9542055   0.9968077   0.9794524
## rep 3        0.9399962   0.9965583   0.9913565   0.9709413   0.9982570   0.9951301
## rep 4        0.9430030   0.9969896   0.9969133   0.9754376   0.9988682   0.9979196

Background or similarity test

The background or similarity test compares the overlap seen between two species’ ENMs to the overlap expected by chance if one or both species was effectively choosing habitat at random from within their broad geographic range. The purpose of this test is to correct for the availability of habitat and ask whether the observed similarity between species or populations is significantly more (or less) than expected given the available set of environments in the regions in which they occur.

NOTE: In order for models to be comparable, both empirical and pseudoreplicate models for the background test are conducted with pseudoabsence points pooled for the two species being compared.

In Warren et al. 2008, we developed this test in the context of comparing one species’ actual occurrence to the random background occurrences of the other species. This is what we call an “asymmetric” test, and in our case we did the test in both directions with the idea that we might compare the results of A vs. B background to the results of B vs. A background. This may be informative in some cases, but many people have also found this asymmetry confusing (and indeed it is often difficult to interpret). For that reason, the background test here can be conducted against a null hypothesis that is generated from “asymmetric” (species.1 vs species.2 background) or “symmetric” (species.1 background vs. species.2 background) comparisons.

Here, for instance, is a Bioclim background test using the classical asymmetric approach:

bg.bc.asym <- background.test(species.1 = monticola, species.2 = cyreni, env = env, type = "bc", nreps = 4, test.type = "asymmetric" )
bg.bc.asym
## 
## 
##  
## 
## Asymmetric background test
##  monticola vs. cyreni background
## 
## background test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.4      0.4      0.4      0.4      0.4      0.4 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  ----------
## empirical    0.0242365   0.1455829   -0.0450249   0.0139750   0.1112424   0.1619493
## rep 1        0.6253512   0.8782090    0.6996700   0.5222534   0.7889548   0.7073589
## rep 2        0.6330398   0.8800637    0.7079561   0.4933189   0.7610055   0.6484262
## rep 3        0.6106309   0.8590911    0.6397875   0.4560589   0.7277260   0.6433115
## rep 4        0.5999843   0.8589928    0.6486787   0.5171324   0.7877419   0.7572882

And here is a Domain background test using the symmetric approach:

bg.dm.sym <- background.test(species.1 = monticola, species.2 = cyreni, env = env, type = "dm", nreps = 4, test.type = "symmetric" )
bg.dm.sym
## 
## 
##  
## 
## Symmetric background test
##  monticola background vs. cyreni background
## 
## background test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.4      0.4      0.4      0.4      0.4      0.4 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  ----------
## empirical    0.1974048   0.4272913   -0.1020590   0.0783117   0.2660741   0.2228439
## rep 1        0.9442736   0.9950074    0.7039206   0.7901624   0.9223426   0.9540107
## rep 2        0.9520233   0.9945641    0.8056707   0.8457025   0.9465923   0.9500261
## rep 3        0.9573928   0.9977192    0.7277525   0.8592384   0.9587948   0.9255293
## rep 4        0.9325482   0.9942441    0.5241340   0.7932810   0.9237870   0.9380176

Ecospat tests

Using enmtools.species objects also provides a simplified interface to the niche equivalency and similarity tests (or identity and background tests, respectively) that were developed by Broennimann et al. (2012). These tests do not rely on ENMs, instead using kernel density smoothing to estimate density of the species in environment space. Ecospat also uses the density of the available environment to correct for availability when measuring overlaps, so that overlaps are not strictly driven by availability of combinations of environmental variables.

These tests only work with two environmental axes, so they are often done with the top two PC axes of a set of environments. In our case we’ll just pick a couple of environmental layers, though (bio1 and bio2). Here’s an equivalency/identity test:

esp.id <- enmtools.ecospat.id(monticola, cyreni, env[[c("bio1", "bio12")]])
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
esp.id
## 
## 
##  
## 
## Ecospat identity test monticola vs. cyreni
## 
## ecospat.id test empirical overlaps:
## $D
## [1] 0.0477392
## 
## $I
## [1] 0.2097866
## 
## 
## 
## ecospat.id test p-values:
## D I 
## 0 0

## NULL

And here’s a symmetric background test. The difference between symmetric and asymmetric for these tests is the same as for the background tests presented above.

esp.bg.sym <- enmtools.ecospat.bg(monticola, cyreni, env[[c("bio1", "bio12")]], test.type = "symmetric")
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
esp.bg.sym
## 
## 
##  
## 
## Ecospat background test symmetric monticola vs. cyreni
## 
## ecospat.bg test empirical overlaps:
## $D
## [1] 0.0477392
## 
## $I
## [1] 0.2097866
## 
## 
## 
## ecospat.bg test p-values:
##    D    I 
## 0.52 0.48

## NULL

Note that if you provide more than two layers to the enmtools.ecospat function, it will performa a PCA analysis on the provided layers and measure overlaps on the first two axes of that PCA space.

esp.bg.sym <- enmtools.ecospat.bg(monticola, cyreni, env, test.type = "symmetric")
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## [1] "More than two layers in environment stack and no layers argument passed, performing PCA..."
esp.bg.sym
## 
## 
##  
## 
## Ecospat background test symmetric monticola vs. cyreni
## 
## ecospat.bg test empirical overlaps:
## $D
## [1] 0.01293428
## 
## $I
## [1] 0.109648
## 
## 
## 
## ecospat.bg test p-values:
##    D    I 
## 0.40 0.34

## NULL

Rangebreak tests

ENMTools also allows you to perform linear, blob, and ribbon rangebreak tests as developed in Glor and Warren 2011. The linear and blob tests are two versions of a test that permit one to ask whether the geographic regions occupied by two species are more environmentally different than expected by chance. The ribbon test, meanwhile, is designed to test whether the ranges of two species are divided by a region that is relatively unsuitable to one or both forms.

For the linear and blob tests, you call them very much like you would the identity and background tests. Here’s a linear one using GLM models:

rbl.glm <- rangebreak.linear(monticola, cyreni, env, type = "glm", nreps = 4)
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## Building empirical models...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Building replicate models...
## 
## Replicate 1 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 2 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 3 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 4 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
rbl.glm
## 
## 
##  
## 
## Linear rangebreak test monticola vs. cyreni
## 
## rangebreak test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.6      0.2      0.4      0.6 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I      env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  -----------
## empirical    0.3244487   0.5929443    0.1049357   0.3693050   0.6167357    0.1377594
## rep 1        0.4809015   0.6995934   -0.0153935   0.4081772   0.5856022   -0.3698372
## rep 2        0.5815106   0.8428696    0.3001375   0.5575885   0.7849459    0.5155767
## rep 3        0.5574340   0.8181974    0.7127636   0.6869428   0.8554630    0.7736962
## rep 4        0.4362774   0.7276420   -0.1933681   0.5250224   0.7514190   -0.0216100

And here’s a blob test using Bioclim:

rbb.bc <- rangebreak.blob(monticola, cyreni, env, type = "bc", nreps = 4)
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## Building empirical models...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Building replicate models...
## 
## Replicate 1 ...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 2 ...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 3 ...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 4 ...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
rbb.bc
## 
## 
##  
## 
## blob rangebreak test monticola vs. cyreni
## 
## rangebreak test p-values:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.2      0.2      0.2      0.2      0.2      0.2 
## 
## 
## Replicates:
## 
## 
##                      D           I     rank.cor       env.D       env.I     env.cor
## ----------  ----------  ----------  -----------  ----------  ----------  ----------
## empirical    0.0242365   0.1455829   -0.0450249   0.0139619   0.1113288   0.1604440
## rep 1        0.2434010   0.4642500    0.3530777   0.1812487   0.3970032   0.2835958
## rep 2        0.4869789   0.7539492    0.7106392   0.3267273   0.5474011   0.4026797
## rep 3        0.7112108   0.9242979    0.8548985   0.5655792   0.8121346   0.7726571
## rep 4        0.4052743   0.5537883    0.3254376   0.2447240   0.3340848   0.3189047

If you want to access the individual replicates (for instance to see how your ranges are being split up), you can find them in the list named “replicate.models” inside your rangebreak test object.

rbl.glm$replicate.models$monticola.rep.1
## 
## 
## Formula:  presence ~ bio1 + bio12 + bio7
## <environment: 0x7fc1a07d9400>
## 
## 
## Data table (top ten lines): 
## 
##  Longitude   Latitude   bio1   bio12   bio7   presence
## ----------  ---------  -----  ------  -----  ---------
##   0.130000   42.64000     61    1053    238          1
##  -0.020000   42.59000     62    1027    240          1
##  -2.450000   42.46000    126     512    265          1
##  -5.132756   43.49572    133     822    190          1
##  -5.160000   43.47000    133     822    190          1
##  -5.160000   43.47000    133     822    190          1
##  -4.941888   43.35310    128     843    194          1
##  -4.790000   43.29000     90     958    223          1
##  -4.790000   43.29000     90     958    223          1
##  -5.290000   43.47000    130     829    193          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9546  -0.7398  -0.3729   0.6885   1.8557  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 15.9363092  2.5883132   6.157 7.41e-10 ***
## bio1        -0.0555566  0.0074661  -7.441 9.98e-14 ***
## bio12       -0.0012044  0.0007862  -1.532    0.126    
## bio7        -0.0367118  0.0057615  -6.372 1.87e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 720.87  on 901  degrees of freedom
## Residual deviance: 571.45  on 898  degrees of freedom
## AIC: 261.27
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 260 
## n absences     : 642 
## AUC            : 0.7887102 
## cor            : 0.4030759 
## max TPR+TNR at : 0.219305 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 260 
## n absences     : 10000 
## AUC            : 0.4294042 
## cor            : 0.002615941 
## max TPR+TNR at : 0.4021146 
## 
## 
## Proportion of data wittheld for model fitting:  0
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 0.004360313, 0.9998061  (min, max)
## 
## 
## 
## Notes:

rbl.glm$replicate.models$cyreni.rep.1
## 
## 
## Formula:  presence ~ bio1 + bio12 + bio7
## <environment: 0x7fc19ff94838>
## 
## 
## Data table (top ten lines): 
## 
##  Longitude   Latitude   bio1   bio12   bio7   presence
## ----------  ---------  -----  ------  -----  ---------
##  -4.010000    40.7800     91     555    284          1
##  -3.860000    40.7200     94     568    284          1
##  -7.010000    41.9400    110     931    265          1
##  -7.010000    41.9400    110     931    265          1
##  -3.890000    40.6900     94     568    284          1
##  -3.890000    40.6900     94     568    284          1
##  -3.890000    40.6900     94     568    284          1
##  -4.014862    40.7364     91     555    284          1
##  -4.130000    40.7800     91     555    284          1
##  -4.130000    40.7800     91     555    284          1
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.79118  -0.35179  -0.22900  -0.03461   1.75918  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -51.520936  10.866974  -4.741 2.13e-06 ***
## bio1         -0.035119   0.010757  -3.265   0.0011 ** 
## bio12         0.011270   0.002483   4.539 5.65e-06 ***
## bio7          0.170618   0.034134   4.998 5.78e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 210.72  on 847  degrees of freedom
## Residual deviance: 139.24  on 844  degrees of freedom
## AIC: 74.252
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 76 
## n absences     : 772 
## AUC            : 0.8574107 
## cor            : 0.2391515 
## max TPR+TNR at : -0.1355184 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 76 
## n absences     : 10000 
## AUC            : 0.7046934 
## cor            : 0.07078544 
## max TPR+TNR at : 0.2127089 
## 
## 
## Proportion of data wittheld for model fitting:  0
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 6.598746e-13, 0.9999969  (min, max)
## 
## 
## 
## Notes:

For the ribbon rangebreak test, you will need one extra thing; a third enmtools.species object representing the occurrence points (for one or both species) that fall within the ribbon of putatively unsuitable habitat. In the case of these two lizards we don’t have such a ribbon, so we’ll just simulate one based on some random points.

plot(env[[1]])
points(cyreni$presence.points, col = "red")
points(monticola$presence.points, col = "blue")

ribbon <- enmtools.species(species.name = "ribbon")
ribbon$presence.points <- data.frame(Longitude = runif(n = 10, min = -9, max = 0),
                                      Latitude = runif(n = 10, min = 40.5, max = 42))
points(ribbon$presence.points, pch = 16)

ribbon$range <- background.raster.buffer(ribbon$presence.points, 20000, mask = env)
ribbon
## 
## 
## Range raster: 
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : bio1 
## values     : 1, 1  (min, max)
## 
## 
## 
## Presence points (first ten only): 
## 
##   Longitude   Latitude
## -----------  ---------
##  -1.4820608   40.63507
##  -8.5952518   41.73363
##  -0.7534716   40.66044
##  -8.8531698   41.02226
##  -2.1542777   41.48831
##  -0.0117659   41.42271
##  -0.3322167   41.74504
##  -7.1280546   41.51917
##  -0.8396206   41.20462
##  -4.9374491   40.96968
## 
## 
## Background points not defined.
## 
## Species name:  ribbon

Now we’ll run a ribbon rangebreak test using GLM models with quadratic effects. We also need to tell it the width of the ribbons to generate for the replicates. The units for the width argument are the same units that the presence points are in; e.g., if the points are in decimal degrees you should supply the width of the barrier in decimal degrees.

rbr.glm <- rangebreak.ribbon(monticola, cyreni, ribbon, env, type = "glm", f = pres ~ poly(bio1, 2) + poly(bio12, 2) + poly(bio7, 2), width = 0.5, nreps = 4)
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## 
## Building empirical models...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Building replicate models...
## 
## Replicate 1 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 2 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 3 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## 
## Replicate 4 ...
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species ribbon 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## Drawing background from species background points.
## 
## Adding environmental data to species outside 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
## [1] "Trying to find starting conditions, attempt 1"
## [1] "Building replicates..."
rbr.glm
## 
## 
##  
## 
## ribbon rangebreak test monticola vs. cyreni
## 
## rangebreak test p-values...
## 
## Species 1 vs. Species 2:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.8      0.4      0.8      0.4      0.4      0.8 
## 
## Species 1 vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      1.2      1.2      0.4      1.2      1.2      0.4 
## 
## Species 2 vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      0.4      0.4      0.4      0.4      0.4      0.4 
## 
## Outside vs. Ribbon:
##        D        I rank.cor    env.D    env.I  env.cor 
##      1.2      1.2      0.4      1.2      1.2      0.4 
## 
## 
## Replicates:
## 
## Species 1 vs. Species 2:
##                    D         I    rank.cor      env.D     env.I      env.cor
## empirical 0.09958614 0.2836189 -0.13381525 0.06046937 0.1532053  0.009763588
## rep 1     0.20363542 0.4435034 -0.21035675 0.23869726 0.3965820 -0.043112346
## rep 2     0.37226018 0.6160198 -0.21315142 0.24355023 0.4227148 -0.112987332
## rep 3     0.37598089 0.5933995 -0.34135732 0.18061195 0.3011252 -0.299113827
## rep 4     0.09691352 0.2910089  0.07892328 0.11998350 0.3115225  0.317139800
## 
## Species 1 vs. Ribbon:
##                    D         I    rank.cor      env.D     env.I     env.cor
## empirical 0.49081076 0.7383465 -0.18782530 0.24585701 0.4471413 -0.05452458
## rep 1     0.09423682 0.2666710 -0.16218057 0.20515732 0.3125827  0.07096497
## rep 2     0.55452653 0.7905449  0.09304785 0.31658589 0.5380501  0.06863749
## rep 3     0.53623312 0.7996006  0.45056886 0.33498857 0.5426246  0.30224947
## rep 4     0.14260851 0.3624106  0.35949853 0.06215738 0.2280846  0.38314114
## 
## Species 2 vs. Ribbon:
##                   D         I   rank.cor      env.D     env.I    env.cor
## empirical 0.1271460 0.2847418 -0.1720820 0.04043136 0.1080127 -0.5676745
## rep 1     0.5322678 0.7860151  0.7842682 0.67232139 0.8413233  0.8656964
## rep 2     0.7689889 0.9527502  0.9112095 0.90997509 0.9769810  0.9739647
## rep 3     0.6501376 0.8226096  0.5791852 0.53248151 0.7246371  0.7784277
## rep 4     0.2954128 0.4951645  0.8364279 0.28070106 0.5123339  0.6146980
## 
## Outside vs. Ribbon:
##                   D         I   rank.cor      env.D     env.I     env.cor
## empirical 0.4753460 0.7184761 -0.4439931 0.26050847 0.4619048 -0.13205995
## rep 1     0.1017739 0.2951013 -0.1555732 0.25713760 0.3989164 -0.00551546
## rep 2     0.6621902 0.8722107  0.3147381 0.44905364 0.6526806  0.13617149
## rep 3     0.6036415 0.8491942  0.6483865 0.45283812 0.6841596  0.54005250
## rep 4     0.1474755 0.3723265  0.3942329 0.05201737 0.2047028  0.30523560

Note that the output table here has slope, intercept, and intercept offset.

rbr.glm$lines.df
##         slope  intercept     offset
## 1   0.4809375   42.87959  0.2774100
## 2 -51.1111363 -207.75568 12.7802295
## 3  -0.7763185   38.32776  0.3164915
## 4   1.1674845   51.20820  0.3843029

The intercept denotes the intercept corresponding to the CENTER of each ribbon. To get the lines denoting the edges of the ribbons (for example if you want to plot the ribbons on a map), you add and substract the offset. In other words, the top edge of the ribbon is given by y = (slope * x) + intercept + offset, while the bottom edge is given by y = (slope * x) + intercept - offset.

Building an enmtools.clade object

Some of the tests in ENMTools, including some really neat ones that are still in development, require you to build an enmtools.clade object. These objects are simply lists that contain a phylogeny and a set of enmtools.species objects. It’s important that the names of the species objects and their species.name attributes match the names in the phylogeny’s tip.labels. For demonstration, we’re going to build an object for a clade of five anoles from Hispaniola. We have the tree, so we’re just going to grab occurrence data from GBIF using the rgbif package.

library(rgbif)
library(ape)

tree.path <- paste(system.file(package="ENMTools"), "/StarBEAST_MCC.species.txt", sep='')

hisp.anoles <- read.nexus(file = tree.path)

keepers <- c("brevirostris", "marron", "caudalis", "websteri", "distichus")

hisp.anoles <- drop.tip(phy = hisp.anoles, tip = hisp.anoles$tip.label[!hisp.anoles$tip.label %in% keepers])
plot(hisp.anoles)

So there’s our tree. Now we’re going to grab some environmental data.

hisp.env <- raster::getData('worldclim', var='bio', res=10)
hisp.env <- raster::crop(hisp.env, extent(-75, -65, 16, 21))

hisp.env <- setMinMax(hisp.env)

And then we’ll create a function to build species from GBIF.

# Automate the process of downloading data and removing duds and dupes
species.from.gbif <- function(genus, species, name = NA, env){

  # Name it after the species epithet unless told otherwise
  if(is.na(name)){
    name <- species
  }

  # Get GBIF data
  this.sp <- enmtools.species(presence.points = gbif(genus = genus, species = species)[,c("lon", "lat")],
                              species.name = name)

  # Rename columns, get rid of duds
  colnames(this.sp$presence.points) <- c("Longitude", "Latitude")
  this.sp$presence.points <- this.sp$presence.points[complete.cases(extract(env, this.sp$presence.points)),]
  this.sp$presence.points <- this.sp$presence.points[!duplicated(this.sp$presence.points),]

  this.sp$range <- background.raster.buffer(this.sp$presence.points, 50000, mask = hisp.env)

  return(this.sp)
}

Now we’ll create five species and add them to a species.clade object that is called brev.clade.

brevirostris <- species.from.gbif(genus = "Anolis", species = "brevirostris", env = hisp.env)
marron <- species.from.gbif(genus = "Anolis", species = "marron", env = hisp.env)
caudalis <- species.from.gbif(genus = "Anolis", species = "caudalis", env = hisp.env)
websteri <- species.from.gbif(genus = "Anolis", species = "websteri", env = hisp.env)
distichus <- species.from.gbif(genus = "Anolis", species = "distichus", env = hisp.env)


brev.clade <- enmtools.clade(species = list(brevirostris, marron, caudalis, websteri, distichus), tree = hisp.anoles)
check.clade(brev.clade)
## 
## 
## An enmtools.clade object with 5 species
## 
## Species names: 
##   brevirostris    caudalis    distichus   marron      websteri
## 
## Tree: 
## 
## Phylogenetic tree with 5 tips and 4 internal nodes.
## 
## Tip labels:
## [1] "brevirostris" "caudalis"     "distichus"    "marron"       "websteri"    
## 
## Rooted; includes branch lengths.
## 
## 
## Data Summary: 
## 
## 
##                species.names   in.tree   presence   background   range   
## -------------  --------------  --------  ---------  -----------  --------
## brevirostris   brevirostris    TRUE      200        0            present 
## caudalis       caudalis        TRUE      26         0            present 
## distichus      distichus       TRUE      834        0            present 
## marron         marron          TRUE      14         0            present 
## websteri       websteri        TRUE      21         0            present

That’s one way to build a clade object by hand, but there’s already one built into ENMTools to experiment with so we’ll just use that.

data(iberolacerta.clade)

Age-overlap correlation tests (AOC)

The AOC tests allow you to examine patterns of range, point, and ENM overlap in the context of a phylogeny. This is effectively a generalized version of several analyses: age-range correlation (e.g., Fitzpatrick and Turelli 2006), ENM overlap in the context of a phylogeny (e.g., Knouft et al. 2006, Warren et al. 2008), and point overlaps (e.g., Cardillo and Warren 2016).

These tests require the creation of an enmtools.clade object, as above. AOC tests consist of two steps: first, the average overlap at each node in the phylogeny is calcualted using a method that takes tree topology into account (see Fitzpatrick and Turelli 2006), then we perform a linear regression to measure the relationship between node age and average overlap. Due to the fact that these overlaps violate many of the assumptions of a regular linear regression, however (e.g., errors are not iid), we can’t calculate significance in the typical way. Instead we performa Monte Carlo test, permuting the identity of the tips of the tree and repeating the node averaging and modeling steps. Finally we measure statistical significance by comparing the empirical slope and intercept to the distribution of slopes and intercepts from the Monte Carlo replicates.

First, let’s do one using geog.range.overlaps, as in Fitzpatrick and Turelli 2006. Note that this analysis requires that each of your species have a range raster stored in their species object (we did that as part of the function used above).

range.aoc <- enmtools.aoc(clade = iberolacerta.clade,  nreps = 50, overlap.source = "range")
summary(range.aoc)
## 
## 
## Age-Overlap Correlation test
## 
## 50 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##       0.03921569       0.03921569

## NULL

Now we can do one using point overlaps just by changing the overlap.source argument:

point.aoc <- enmtools.aoc(clade = iberolacerta.clade,  nreps = 50, overlap.source = "points")
summary(point.aoc)
## 
## 
## Age-Overlap Correlation test
## 
## 50 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##        0.2745098        0.4313725

## NULL

Or we can use similarity between ENMs built for each species. Here we’ll use GLM models:

glm.aoc <- enmtools.aoc(clade = iberolacerta.clade,  env = env, nreps = 50, overlap.source = "glm", f = pres ~ poly(bio1, 2) + poly(bio12, 2))
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species monticola 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species martinezricai 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species cyreni 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species horvathi 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species aurelioi 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species aranica 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
## 
## 
## 
## No background points provided, drawing background from range raster.
## 
## Adding environmental data to species bonnali 
##  Processing presence points...
##  Processing background points...
## 
## 
## Drawing background from species background points.
summary(glm.aoc)
## 
## 
## Age-Overlap Correlation test
## 
## 50 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##        0.2352941        0.3137255

## NULL

Literature cited

Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz, N. G., Thuiller, W., Fortin, M.-J., Randin, C., Zimmermann, N. E., Graham, C. H. and Guisan, A. (2012), Measuring ecological niche overlap from occurrence and spatial environmental data. Global Ecology and Biogeography, 21: 481–497. doi:10.1111/j.1466-8238.2011.00698.x

Fitzpatrick, B. M., & Turelli, M. (2006). The geography of mammalian speciation: mixed signals from phylogenies and range maps. Evolution, 60(3), 601-615.

Knouft, J. H., Losos, J. B., Glor, R. E., & Kolbe, J. J. (2006). Phylogenetic analysis of the evolution of the niche in lizards of the Anolis sagrei group. Ecology, 87(sp7).

Levins, R. 1968. Evolution In Changing Environments. Monographs in Population Biology, volume 2. Princeton University Press, Princeton, New Jersey, USA.

Schoener, T. W. 1968. Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49:704- 726.

Warren, D.L., R.E. Glor, and M. Turelli. 2008. Environmental niche identity versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868-2883. doi: 10.1111/j.1558-5646.2008.00482.x

Warren, D.L., M. Cardillo, D.F. Rosauer, and D.I. Bolnick. 2014. Mistaking geography for biology: inferring processes from species distributions. Trends in Ecology and Evolution 29 (10), 572-580. doi: 10.1016/j.tree.2014.08.003

Warren, D.L., L. Beaumont, R. Dinnage, and J. Baumgartner. 2019. New methods for measuring ENM breadth and overlap in environmental space. Ecography. Doi: 10.1111/ecog.03900